SOA Research Institute Just In Time AWE - LA County Wildfires 202501
https://www.climatologylab.org/gridmet.html
Link to the data files
https://www.northwestknowledge.net/metdata/data/vpd_2025.nc
https://www.northwestknowledge.net/metdata/data/vpd_2024.nc
https://www.northwestknowledge.net/metdata/data/vpd_2023.nc
https://www.northwestknowledge.net/metdata/data/vpd_2022.nc
https://www.northwestknowledge.net/metdata/data/vpd_2021.nc
https://www.northwestknowledge.net/metdata/data/vpd_2020.nc
https://www.northwestknowledge.net/metdata/data/vs_2025.nc
https://www.northwestknowledge.net/metdata/data/bi_2025.nc
https://www.northwestknowledge.net/metdata/data/bi_2024.nc
InĀ [2]:
#!pip install dask
#!pip install xarray
import xarray as xr
import dask
import matplotlib.pyplot as plt
from netCDF4 import Dataset as NetCDFFile
import numpy as np
from mpl_toolkits.basemap import Basemap
InĀ [3]:
vpd = xr.open_mfdataset('vpd_*.nc',combine = 'by_coords')
vpd
Out[3]:
<xarray.Dataset> Size: 12GB
Dimensions: (day: 1841, lat: 585, lon: 1386, crs: 1)
Coordinates:
* lon (lon) float64 11kB -124.8 -124.7 ... -67.06
* lat (lat) float64 5kB 49.4 49.36 ... 25.11 25.07
* day (day) datetime64[ns] 15kB 2020-01-01 ... 202...
* crs (crs) uint16 2B 3
Data variables:
mean_vapor_pressure_deficit (day, lat, lon) float64 12GB dask.array<chunksize=(61, 98, 231), meta=np.ndarray>
Attributes: (12/22)
geospatial_bounds_crs: EPSG:4326
Conventions: CF-1.6
geospatial_bounds: POLYGON((-124.7666666333333 49.40000000000000...
geospatial_lat_min: 25.066666666666666
geospatial_lat_max: 49.40000000000000
geospatial_lon_min: -124.7666666333333
... ...
last_permanent_slice: 306
last_early_slice: 366
last_provisional_slice: 360
note3: Data in slices after last_permanent_slice (1-...
note4: Data in slices after last_provisional_slice (...
note5: Days correspond approximately to calendar day...InĀ [4]:
#the average wind speed measured at a height of 10 meters above the ground, calculated over a 24-hour period
wind_speed = xr.open_mfdataset('vs_*.nc',combine = 'by_coords')
wind_speed
Out[4]:
<xarray.Dataset> Size: 2GB
Dimensions: (day: 383, lat: 585, lon: 1386, crs: 1)
Coordinates:
* lon (lon) float64 11kB -124.8 -124.7 -124.7 ... -67.14 -67.1 -67.06
* lat (lat) float64 5kB 49.4 49.36 49.32 49.28 ... 25.15 25.11 25.07
* day (day) datetime64[ns] 3kB 2024-01-01 2024-01-02 ... 2025-01-17
* crs (crs) uint16 2B 3
Data variables:
wind_speed (day, lat, lon) float64 2GB dask.array<chunksize=(61, 98, 231), meta=np.ndarray>
Attributes: (12/19)
geospatial_bounds_crs: EPSG:4326
Conventions: CF-1.6
geospatial_bounds: POLYGON((-124.7666666333333 49.40000000000000...
geospatial_lat_min: 25.066666666666666
geospatial_lat_max: 49.40000000000000
geospatial_lon_min: -124.7666666333333
... ...
date: 18 January 2025
note1: The projection information for this file is: ...
note2: Citation: Abatzoglou, J.T., 2013, Development...
note3: Data in slices after last_permanent_slice (1-...
note4: Data in slices after last_provisional_slice (...
note5: Days correspond approximately to calendar day...InĀ [5]:
def filter_data_for_la(dataset, min_lat=33.5, max_lat=34.5, min_lon=-118.5, max_lon=-117.5):
lat = dataset['lat']
lon = dataset['lon']
# Filter the data within the bounding box
filtered_data = dataset.where((lat >= min_lat) & (lat <= max_lat) & (lon >= min_lon) & (lon <= max_lon), drop=True)
return filtered_data
InĀ [6]:
vpd_LA = filter_data_for_la(vpd)
wind_speed_LA = filter_data_for_la(wind_speed)
InĀ [7]:
map = Basemap(projection='merc',llcrnrlon=-121.,llcrnrlat=33.,urcrnrlon=-117.,urcrnrlat=36.,resolution='i')
InĀ [8]:
lat = vpd.sel(day = "2025-01-08")['lat']
lon = vpd.sel(day = "2025-01-08")['lon']
lons,lats= np.meshgrid(lon,lat)
x,y = map(lons,lats)
snapshot1 = map.contourf(x,y,vpd.sel(day = "2025-01-08")['mean_vapor_pressure_deficit'],cmap='RdGy_r')
map.colorbar(snapshot1,"bottom", size="5%", pad="2%").set_label('VPD (Kpa)')
plt.title('2025-01-08 mean vapor pressure deficit', fontsize=10)
map.drawcounties()
# text color = 'magenta'
x1, y1 = map(-118.7, 34.7) # Convert longitude, latitude to map coordinates
plt.text(x1, y1, 'Los Angeles County', fontsize=6, color='Navy',fontweight='bold')
x2, y2 = map(-119, 34.05) # Convert longitude, latitude to map coordinates
plt.text(x2, y2, 'Palisades Fire', fontsize=6, color='Navy',fontweight='bold')
x3, y3 = map(-118.4, 34.2) # Convert longitude, latitude to map coordinates
plt.text(x3, y3, 'Eaton Fire', fontsize=6, color='Navy',fontweight='bold')
Out[8]:
Text(289106.67314017046, 160205.35287836706, 'Eaton Fire')
InĀ [9]:
lat = wind_speed.sel(day = "2025-01-08")['lat']
lon = wind_speed.sel(day = "2025-01-08")['lon']
lons,lats= np.meshgrid(lon,lat)
x,y = map(lons,lats)
snapshot2 = map.contourf(x,y,wind_speed.sel(day = "2025-01-08")['wind_speed'],cmap='RdGy_r')
map.colorbar(snapshot2,"bottom", size="5%", pad="2%").set_label('m/s')
plt.title('2025-01-08 Daily Mean Wind Speed (10 m)',fontsize=10)
map.drawcounties()
x1, y1 = map(-118.7, 34.7) # Convert longitude, latitude to map coordinates
plt.text(x1, y1, 'Los Angeles County', fontsize=6, color='Navy',fontweight='bold')
Out[9]:
Text(255748.2108547664, 227627.6431298973, 'Los Angeles County')
InĀ [10]:
vpd_LA.mean_vapor_pressure_deficit.mean(dim=["lat", "lon"]).plot()
plt.title(f'LA county: mean_vapor_pressure_deficit')
plt.xlabel('Day')
#vpd_LA.mean_vapor_pressure_deficit.max(dim=["lat", "lon"]).plot()
#plt.title(f'LA county: Maximum mean_vapor_pressure_deficit')
#plt.xlabel('Day')
Out[10]:
Text(0.5, 0, 'Day')
InĀ [11]:
import numpy as np
from scipy.ndimage import gaussian_filter1d
def gaussian_smoothing(time_series, sigma=1):
"""
Applies Gaussian smoothing to a raw time series.
Parameters:
time_series (array-like): The raw time series data.
sigma (float): The standard deviation for Gaussian kernel. Default is 1.
Returns:
smoothed_series (array-like): The smoothed time series.
"""
smoothed_series = gaussian_filter1d(time_series, sigma=sigma)
return smoothed_series
vpd_mean = vpd_LA.mean_vapor_pressure_deficit.mean(dim=["lat", "lon"])
data_df = vpd_mean.to_dataframe().reset_index()
raw_time_series = data_df['mean_vapor_pressure_deficit']
smoothed_series = gaussian_smoothing(raw_time_series, sigma=20)
#print("Raw time series:", raw_time_series)
print("Smoothed time series:", smoothed_series)
plt.figure(dpi=300,figsize=(12,3))
plt.plot(data_df['day'],data_df['mean_vapor_pressure_deficit'])
plt.plot(data_df['day'],smoothed_series)
plt.grid()
plt.xlabel('Day')
plt.ylabel('Mean vapor pressure deficit')
plt.title(f'Los Angeles county: mean_vapor_pressure_deficit', fontsize=10)
Smoothed time series: [0.82039447 0.82055062 0.82086192 ... 0.89924112 0.89876858 0.89853465]
Out[11]:
Text(0.5, 1.0, 'Los Angeles county: mean_vapor_pressure_deficit')
Appendix: Not used in the report
InĀ [12]:
# this is not used in the report
burn_index = xr.open_mfdataset('bi_*.nc',combine = 'by_coords')
burn_index
Out[12]:
<xarray.Dataset> Size: 2GB
Dimensions: (day: 383, lat: 585, lon: 1386, crs: 1)
Coordinates:
* lon (lon) float64 11kB -124.8 -124.7 -124.7 ... -67.1 -67.06
* lat (lat) float64 5kB 49.4 49.36 49.32 ... 25.15 25.11 25.07
* day (day) datetime64[ns] 3kB 2024-01-01 ... 2025-01-17
* crs (crs) uint16 2B 3
Data variables:
burning_index_g (day, lat, lon) float64 2GB dask.array<chunksize=(61, 98, 231), meta=np.ndarray>
Attributes: (12/22)
geospatial_bounds_crs: EPSG:4326
Conventions: CF-1.6
geospatial_bounds: POLYGON((-124.7666666333333 49.40000000000000...
geospatial_lat_min: 25.066666666666666
geospatial_lat_max: 49.40000000000000
geospatial_lon_min: -124.7666666333333
... ...
last_permanent_slice: 306
last_early_slice: 366
last_provisional_slice: 360
note3: Data in slices after last_permanent_slice (1-...
note4: Data in slices after last_provisional_slice (...
note5: Days correspond approximately to calendar day...InĀ [13]:
#this is not used in the report
lat = burn_index.sel(day = "2025-01-08")['lat']
lon = burn_index.sel(day = "2025-01-08")['lon']
lons,lats= np.meshgrid(lon,lat)
x,y = map(lons,lats)
snapshot3 = map.contourf(x,y,burn_index.sel(day = "2025-01-08")['burning_index_g'],cmap='RdGy_r')
map.colorbar(snapshot3,"bottom", size="5%", pad="2%").set_label('Unitless')
plt.title('2025-01-08 burning_index_g')
map.drawcounties()
Out[13]:
<matplotlib.collections.PolyCollection at 0x127577a10>
InĀ [14]:
#this is not used in the report
#sporadic in winter, more consistent in summer
burn_index_LA = filter_data_for_la(burn_index)
burn_index_LA.burning_index_g.mean(dim=["lat", "lon"]).plot()
#burn_index_LA.burning_index_g.max(dim=["lat", "lon"]).plot()
Out[14]:
[<matplotlib.lines.Line2D at 0x130be5cd0>]
